d <- read.csv(file="./wvfn.csv",head=FALSE,sep=",")
t <- d$V1; x <- d$V2; rho <- d$V4

library(akima)
ex = seq(min(t),max(t),length=200)
why = seq(min(x),max(x),length=200)
zz <- interp(t,x,-rho,xo=ex,yo=why)


data <- read.csv("./RichardsonRK4.csv",head=FALSE,sep=",")
# now can access each as data$V1 or data$V2


pdf("./OrbitsOnContour.pdf") # makes a pdf printing device for plots

filled.contour(zz,xlab="Time",ylab="Position",ylim=c(-10,10),col=rainbow(20,start=0/6,end=4/6),plot.axes = {points(data$V1,data$V2, pch=".",cex=0.6); axis(1); axis(2) })

dev.off() #closes the pdf printing device